Ansatz from Non-Linear Optics Applied to Trapped Bose-Einstein Condensates 
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A simple analytical ansatz, which has been used to describe the intensity profile of the similariton 
laser P, U is used as a variational wave function to solve the Gross-Pitaevskii equation for a wide 
range of interaction parameters. The variational form interpolates between the non-interacting 
density profile and the strongly interacting Thomas-Fermi profile smoothly. The simple form of the 
ansatz is modified for both cylindrically symmetric and completely anisotropic harmonic traps. The 
resulting ground state density profile and energy are in very good agreement with both the analytical 
solutions in the limiting cases of interaction and the numerical solutions in the intermediate regime. 
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It is common for dynamical systems with weak cou- 
pling to show Gaussian behavior with respect to key 
parameters, whereas strongly coupled and highly non- 
linear systems tend to exhibit power law dependencies. 
A striking recent example of this is the ultrashort pulse 
formation in the similariton laser Another nonlin- 
ear system, which exhibits a parabolic density profile for 
strong interactions is a trapped Bose-Einstein Conden- 
sate (BEC). In Bose-Einstein condensation, the density 
of the condensate is analogous to the intensity of light 
in nonlinear optics and the nonlinear governing equation 
for this system, which is called the Gross-Pitaevskii equa- 
tion (GPE) is very similar to the equation for propaga- 
tion of laser light in a nonlinear optical medium. Based 
on this similarity, it is natural to expect similar solu- 
tions for these vastly different systems. Soliton-like and 
self-similar solutions of nonlinear Schrodinger Equation 
(NLSE) are important both in BEC and NLO. 

Soliton-like solutions arise when the nonlinearity is 
compensated by the dispersion, and they are the only 
exact analytical solutions of these NLSEs, whereas self- 
similar solutions are asymptotic solutions that show up 
when the effects of initial conditions die out but the sys- 
tem is still far from the final state 0. Although soliton 
type solutions have been extensively studied in both the 
NLO and BEC communities, self-similar solutions are not 
as comprehensively investigated. In optics, these type of 
solutions are used more extensively from Raman scatter- 
ing to pulse propagation in fibers, and it is shown that 
linearly chirped parabolic pulses are exact asymptotic so- 
lutions of NLSE with gain Q . Recently, we have demon- 
strated experimentally and numerically that self-similar 
propagation of ultrashort parabolic pulses (similaritons) 
are stable in a laser resonator More recently, we 
have developed a semi-analytic theory of the similariton 
laser 0, Q . Instrumental in this step was the introduc- 
tion of an ansatz to describe the intensity profile of this 
pulse, which can be 'tuned' to any condition, ranging 
from weakly nonlinear (Gaussian pulses) to strongly non- 
linear (parabolic pulses). Here, motivated by the math- 
ematical similarity between the two systems, we apply 
the same ansatz to describe the density profile of a BEC 
in a quadratic trap. We show that this ansatz describes 
the system with excellent accuracy, in the whole range 
from the non-interacting limit to the strongly interacting 



limit. 

Gross-Pitaevskii theory [G] gives a very successful de- 
scription of the ground state and excitations of the BECs 
in dilute atomic gases. The success of this theory implies 
that the condensate can be described accurately with a 
single wave function and the interactions between the 
particles are through s-wave scatterings. The interac- 
tion of the particles are then represented by the interac- 
tion strength g — "^^^ ° , where a is the s-wave scattering 
length, and m is the mass of the trapped particles. The 
theory reduces to a single equation that describes the 
condensate wave function, known as GPE; a type of non- 
linear Schrodinger equation which arises in many areas 
of physics like NLO and hydrodynamic theory of fluids. 



TABLE L The values of the wave function at the center, the 
root mean square si zes rrm s and chemical potentials are tab- 
ulated in units of ^ N/ , and hcu respectively. For com- 
parison numerical results of Ref. [llj are given in parentheses. 
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^(0) 








0.2496 
0.9986 
2.4964 
9.9857 
49.926 
249.64 
2496.4 


0.4238 (0.4238) 
0.3969 (0.3843) 
0.3475 (0.3180) 
0.2515 (0.2581) 
0.1739 (0.1738) 
0.1097 (0.1066) 
0.0665 (0.0655) 
0.0330 (0.0328) 


1.2248 (1.2248) 
1.2794 (1.2785) 
1.3981 (1.3921) 
1.5355 (1.5356) 
1.8822 (1.8821) 
2.5071 (2.5057) 
3.4152 (3.4145) 
5.3855 (5.3852) 


1.5000 (1.5000) 
1.6805 (1.6774) 
2.0885 (2.0650) 
2.5803 (2.5861) 
4.0089 (4.0141) 
7.2576 (7.2484) 
13.559 (13.553) 
33.812 (33.810) 



GPE can be obtained by minimizing the ground state 
energy functional of the condensate 



E{^|;) = / dr 
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with respect to the wave function. The terms in the en- 
ergy functional correspond to kinetic, trapping and inter- 
action energies, respectively. The trapping potential can 
generally be approximated with a harmonic potential for 
many of the experiments. Time-independent GPE fol- 
lows as, 
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V2^(r) -I- V{r)ij{r) + g|^(r)| V(r) = ^i^b{r), (2) 
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FIG. 1; Change of variational parameters with interaction 
parameter (3 is given. Left plot shows the number of terms in 
the summation and the right one shows the change of width 
of the similariton ansatz in units of oscillator length a^j. 
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FIG. 2: (Color online) Wave function calculated with the 
steepest descent method Q is shown with bold line (blue) 
whereas the similariton ansatz solution is given with dashed 
line (red) for (3 = 100. 



where ^ is chemical potential introduced as the Lagrange 
multiplier for the normalization constraint J c?rn(r) = N , 
where n(r) = |?/'(r)P is the density of the condensate. 

Nonlinearity of GPE is due to interaction between par- 
ticles and its effect becomes more pronounced as the 
number of particles in the condensate increases which 
is the case for current experiments where more than 10^ 
particles form the BEC. Since very few exact solutions 
of NLSE, such as solitons, are known , many numeri- 
cal algorithms 0, H_ 1, [l^, [HI and variational methods 
[l3,[l3|,[l3,[l^[ll,[T7| are developed to find ground state 
solutions. Although variational methods give only an up- 
per bound to the ground state energy, they require less 
calculation and can give accurate results if a suitable trial 
function is chosen. Another advantage of the variational 
principle is that it gives the functional form of the wave 
function which can be used to obtain other properties of 
the condensate. 



FIG. 3: (Color online) Ground state energy with respect to in- 
teraction parameter /3 obtained with the variational function 
(black solid line). The resulting energy of a Gaussian vari- 
ational function is given with dotted (red) line and energy 
obtained with Thomas-Fermi solution is given with dashed 
line (green). Improved Thomas- Fermi solution [l^ is given 
with dotted line (blue). Energy per particle is given in units 
of fto). The inset is given for small /? values. 



Therefore, many trial functions are proposed to ob- 
tain a lower bound for the ground state energy. Trial 
functions are generally chosen by adding parameters to a 
known approximate analytical solution. The trivial ap- 
proximate solutions are obtained by looking at the lim- 
iting cases of the GPE where the nonlinearity is negli- 
gibly small or very high. Assuming an ideal Bose gas 
equation ^ reduces to the Schrodinger equation where 
the chemical potential corresponds to the energy eigen- 
value, by neglecting the nonlinear term. With a harmonic 
trap the problem turns into a harmonic oscillator prob- 
lem where the solution is a Gaussian. For the opposite 
case where nonlinearity is dominant, the kinetic energy 
term can be neglected in GPE, and then the equation 
can easily be solved for the density,n(r) = ^~2g^^ ,when 
the right hand side of the equation is positive. This ap- 
proximation is known as Thomas-Fermi approximation 
(TEA) and it shows that as the interaction increases, the 
density profile changes from a Gaussian to a parabola 
for a harmonic trap. TEA can be improved by adding 
the kinetic energy term obtained with the resulting wave 
function with a suitable cut off nm as, 



^ =(15/3)-2/5(i In (480/3)-^). (3) 



With this insight many trial functions are proposed, 
for isotropic traps and anisotropic traps with cylindrical 
symmetry to describe the intermediate regime where nei- 
ther TEA and nor the ideal gas approximation is valid. 
Here, we make use of a simple analytic function which 
has already been used in NLO where a similar behavior 
-Gaussian to parabolic- for intensity profile of a similari- 
ton laser. 

The so-called 'similariton' pulse in optics has a nearly 
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parabolic intensity profile to reduce the effect of Kerr 
nonlinearity. However, if the nonlinearity of the system 
is reduced, the pulse assumes the well-known Gaussian 
shape of the dispersion-managed solitons. Therefore, the 
ansatz proposed in Q to describe these pulses has an 
adjustable profile between a Gaussian and an inverted 
parabola 



Sn{x) = exp (- 



k=l 



This function becomes a Gaussian when ri = 1 



5'„=i(a;) = exp(-a;^), 



(4) 



(5) 



and turns into an inverted parabola when n goes to oo 
since the summation in the exponent converges to ln(l — 



for \x\ < 1, 
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Sn^oo = cxp (- ^ ^) = exp [ln(l 
fe=i 



(6) 

Moreover the expansion converges so quickly that adding 
about ten terms is enough to get a parabolic profile with 
smooth ends. Besides, this function is easily integrable 
which makes it a good candidate for variational calcula- 
tions. 

Motivated by these properties, we use the similariton 
ansatz Eq. |4] as our trial wave function to minimize 
the energy functional given in Eq. (H]). To simplify the 
calculations, we non-dimensionalize the Gross-Pitaevskii 
functional by scaling length, energy and wave function 

with oscillator length a^^ — ^ Hlu and \/Na^, re- 
spectively. We first analyze the solution for spherical 



harmonic trap V{r) = 



and introduce the pa- 



rameter (3 = Na/a^ which is a measure of the strength 
of the interaction. With this rescaling the energy func- 
tional becomes 



N 



2tt / d^r (IW'WP + T/(r)|V'(r)|2 +27r/3|i/;(r)|4) 

.(7) 

Ideally (3 can take any value between —oo to oo since 
all the parameters are experimentally tunable. However, 
negative scattering length which means attractive inter- 
action, causes collapse of the condensate when the parti- 
cle number is high. In this regime, our results agree with 
In the present work, we concentrate on repulsive 
interaction. With proper normalization the trial wave 
fimction has the form, 



1 ^ i^r/dr 

k=l 



(8) 



where d and n are our variational parameters with / = 
Jq°° drr^Snir) which is an integral that can be calculated 
analytically for n = 1, 2 and numerically for n > 2. Here 



the parameter d is responsible for the width of the con- 
densate which increases as the interaction increases, and 
n takes care of flattening of the central density. We min- 
imize the energy with respect to d for different n values 
and chose the n that gives the minimum energy. For d, 
we obtain a fifth order polynomial equation where only 
one of the roots is physically meaningful. 

We compare our results with the analytical approxi- 
mations as follows. For small (3 values, our trial func- 
tion reduces to a Gaussian and gives the exact result for 
/3 = 0, and for large /3 our results agree well with the 
improved TEA results as shown in Eig. [3l We also com- 
pare the resulting wave function with the numerical so- 
lutions obtained by steepest descent method for different 
13 values in Eig. [2l We also tabulate our results in and 
include the results of a recent numerical analysis which 
minimizes the energy functional directly by the finite el- 
ement method. Here it should be noted that tabulated 
kinetic, trap, and interaction energies satisfy the virial 
theorem 2Ekin + ^Epot — 'iEi„t = 0, to our numerical ac- 
curacy. It is also remarkable that even for large /3 adding 
10 terms is enough to find a good approximation for the 
wave function (see Eig. [T]) which shows the simplicity of 
the calculations. 

Using similar trial functions, we can also solve the GPE 
for a cylindrical trap and a fully anisotropic trap. Eor 
the cylindrically symmetric trap, trial function takes the 
form, 



z) = C exp(- 



(p/dp 



l2fc 



fc=l 



2k 



I exp(- 



jz/d 

k=l 



\2k 



2k 



(9) 

^here, C = .f^^^^^Jp - pdpSnip) and I, = 

dzSn{z). We have four variational parameters, but 
calculations are similar to the isotropic case. We compare 
our results with the numerical results of Dalfovo et al. 
in Table 2. Cylindrically symmetric traps are the most 
common traps in BEC setups and aspect ratio obtained 
from is very important to identify the BEC phase 
in these experiments. It is shown in that for the 

noninteracting case this ratio is equal to a/A, and goes to 
A in the Thomas-Eermi limit. This result is clearly seen 
from the values in Table 2 where A — VS and it is also 
evident that convergence of TEA is very slow. 

There are also experiments with fully anisotropic traps 
and for this case the trial function can also be mod- 
ified similarly with six variation parameters. The results 
are given in Table 3, where agreement with the results in 
[Tot is apparent. 

In summary, an ansatz is introduced to investigate the 
ground state properties of a BEC at zero temperature for 
quadratic traps with arbitrary anisotropy. The ground 
state energy and wave function are found to be very 
accurate anywhere from the noninteracting case to the 
highly repulsive one, as compared with numerical stud- 
ies. The form of the trial function changes from a Gaus- 
sian to a parabola smoothly, and it successfully describes 
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TABLE II: Results ol our calculation for a cylindrically harmonic trap with A — \/8- Energy and length units are Ntuj and 
a^,. The results of the numerical calculation in Ref. [0] are given in parentheses for comparison except for the last row. For 
/3 — 2165 Thomas-Fermi result for the chemical potential is given in parentheses. 



p 


Xrms 


Zrms 




Etr 


Eint 


M 


0.0000 


0.7071 (0.707) 


0.4204 (0.42) 


1.2071 (1.207) 


1.2071 (1.207) 


0.0000 (0.000) 


2.4142 


(2.414) 


0.4330 


0.7901 (0.79) 


0.4374 (0.44) 


1.0539 (1.06) 


1.3894 (1.39) 


0.2237 (0.21) 


2.8907 


(2.88) 


0.8660 


0.8500 (0.85) 


0.4472 (0.45) 


0.9976 (0.98) 


1.5225 (1.52) 


0.3500 (0.36) 


3.2200 


(3.21) 


2.1650 


0.9657 (0.96) 


0.4707 (0.47) 


0.8528 (0.86) 


1.8188 (1.81) 


0.6440 (0.63) 


3.9596 


(3.94) 


4.3300 


1.0892 (1.08) 


0.4966 (0.50) 


0.7337 (0.76) 


2.1730 (2.15) 


0.9595 (0.96) 


4.8258 


(4.77) 


8.6600 


1.2319 (1.23) 


0.5332 (0.53) 


0.6709 (0.66) 


2.6549 (2.64) 


1.3227 (1.32) 


5.9712 


(5.93) 


21.650 


1.4798 (1.47) 


0.5930 (0.59) 


0.5314 (0.54) 


3.5963 (3.57) 


2.0432 (2.02) 


8.2142 


(8.14) 


43.300 


1.7038 (1.69) 


0.6536 (0.65) 


0.4351 (0.45) 


4.6121 (4.57) 


2.7847 (2.74) 


10.616 


(10.5) 


64.950 


1.8447 (1.84) 


0.6989 (0.70) 


0.4128 (0.41) 


5.3569 (5.31) 


3.2960 (3.26) 


12.361 


(12.2) 


86.600 


1.9562 (1.94) 


0.7319 (0.73) 


0.3789 (0.38) 


5.9693 (5.91) 


3.7270 (3.68) 


13.802 


(13.7) 


2165 


3.7367 


1.3297 


0.1459 


21.035 


13.926 


49.033 


(48.329) 



TABLE III: The chemical potentials per particle in units of 
fwj are calculated using an ansatz modified for a completely 
anisotropic trap with A — \/2 and (77 = 2. The values of the 
interaction parameter j3 are obtained from Ref. [l^ and the 
values in parentheses corresponds to the numerical solution 
given in [10|. The variational parameters nx,y,z and dx,y,z are 
given. 



p 


n 


X 7 dx 


Wy, dy 


nz,dz 







1 , 


1.000 


1 , 0.840 


1 , 0.707 


2.207 (2.207) 


1.787 


2 , 


1.753 


2 , 1.368 


1 , 0.817 


3.604 (3.572) 


3.575 


2 , 


1.956 


2 , 1.489 


2 , 1.160 


4.385 (4.345) 


7.151 


3 , 


2.433 


2 ,1.654 


2 , 1.258 


5.492 (5.425) 


14.302 


3 , 


2.780 


3 , 2.035 


2 , 1.384 


7.010 (6.904) 


28.605 


4 , 


3.310 


3 , 2.30 


3, 1.687 


9.049 (8.900) 


57.211 


5 , 


3.880 


4 , 2.725 


3 , 1.896 


11.78 (11.57) 



the intermediate regime of moderate interaction. Impor- 
tant quantities like aspect ratio, chemical potential, and 



root mean square size of the clouds are calculated and 
compared to numerical studies [l^, [HI ■ With a sHght 
modification the suggested form of the wave function can 
be applied to vortex states of the condensate. Time- 
dependent GPE can also be solved using a similar form 
to study growth dynamics. 

We have shown previously that a nonhnear system in 
optics, namely a high-energy femtosecond laser oscillator, 
exists stably between two extreme hmits, correspond- 
ing to Gaussian pulse profiles for weak nonlinearity and 
parabohc profiles for strong nonlinearity. This indeed 
appears to be a common behavior observed in many sys- 
tems, including the trapped BEG analyzed in this paper. 
We have a simple analytical function, which has the cru- 
cial property of interpolating any state between these ex- 
tremes. In this paper, we have shown that the variational 
approach with the same ansatz yields excellent results for 
the BEG system in a quadratic trap. We believe that our 
approach can be generally applicable to other nonlinear 
systems in disparate fields. 
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